home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / PIBSIGS / GAMMAIN.PAS < prev    next >
Pascal/Delphi Source File  |  1986-06-22  |  6KB  |  189 lines
  1. (*-------------------------------------------------------------------------*)
  2. (*                GammaIn -- Incomplete Gamma integral                     *)
  3. (*-------------------------------------------------------------------------*)
  4.  
  5. FUNCTION GammaIn(     Y, P    : REAL;
  6.                       Dprec   : INTEGER;
  7.                       MaxIter : INTEGER;
  8.                   VAR Cprec   : REAL;
  9.                   VAR Iter    : INTEGER;
  10.                   VAR Ifault  : INTEGER ) : REAL;
  11.  
  12. (*-------------------------------------------------------------------------*)
  13. (*                                                                         *)
  14. (*       Function:  GammaIn                                                *)
  15. (*                                                                         *)
  16. (*       Purpose:   Evaluates Incomplete Gamma integral                    *)
  17. (*                                                                         *)
  18. (*       Calling Sequence:                                                 *)
  19. (*                                                                         *)
  20. (*            P     := GammaIn(       Y, P    : REAL;                      *)
  21. (*                                    Dprec   : INTEGER;                   *)
  22. (*                                    MaxIter : INTEGER;                   *)
  23. (*                               VAR  Cprec   : REAL;                      *)
  24. (*                               VAR  Iter    : INTEGER;                   *)
  25. (*                               VAR  Ifault  : INTEGER ) : REAL;          *)
  26. (*                                                                         *)
  27. (*                 Y      --- Gamma distrib. value                         *)
  28. (*                 P      --- Degrees of freedom                           *)
  29. (*                 Ifault --- error indicator                              *)
  30. (*                                                                         *)
  31. (*                 P      --- Resultant probability                        *)
  32. (*                                                                         *)
  33. (*       Calls:                                                            *)
  34. (*                                                                         *)
  35. (*            None                                                         *)
  36. (*                                                                         *)
  37. (*       Remarks:                                                          *)
  38. (*                                                                         *)
  39. (*            Either an infinite series summation or a continued fraction  *)
  40. (*            approximation is used, depending upon the argument range.    *)
  41. (*                                                                         *)
  42. (*-------------------------------------------------------------------------*)
  43.  
  44. CONST
  45.    Oflo    = 1.0E+37;
  46.    MinExp  = -87.0;
  47.  
  48. VAR
  49.    F:        REAL;
  50.    C:        REAL;
  51.    A:        REAL;
  52.    B:        REAL;
  53.    Term:     REAL;
  54.    Pn:       ARRAY[1..6] OF REAL;
  55.    Gin:      REAL;
  56.    An:       REAL;
  57.    Rn:       REAL;
  58.    Dif:      REAL;
  59.    Eps:      REAL;
  60.    Done:     BOOLEAN;
  61.  
  62. LABEL 9000;
  63.  
  64. BEGIN (* GammaIn *)
  65.                                    (* Check arguments *)
  66.    Ifault  := 1;
  67.    GammaIn := 1.0;
  68.  
  69.    IF( ( Y <= 0.0 ) OR ( P <= 0.0 ) ) THEN GOTO 9000;
  70.  
  71.                                    (* Check value of F *)
  72.    Ifault := 0;
  73.    F      := P * LN( Y ) - ALGama( P + 1.0 ) - Y;
  74.    IF ( F < MinExp ) THEN GOTO 9000;
  75.  
  76.    F := EXP( F );
  77.    IF( F = 0.0 ) THEN GOTO 9000;
  78.  
  79.                                   (* Set precision *)
  80.    IF Dprec > MaxPrec THEN
  81.       Dprec := MaxPrec
  82.    ELSE IF Dprec <= 0 THEN
  83.       Dprec := 1;
  84.  
  85.    Cprec  := Dprec;
  86.  
  87.    Eps    := PowTen( -Dprec );
  88.  
  89.                                    (* Choose infinite series or *)
  90.                                    (* continued fraction.       *)
  91.  
  92.    IF( ( Y > 1.0 ) AND ( Y >= P ) ) THEN
  93.  
  94.       BEGIN (* Continued Fraction *)
  95.  
  96.          A       := 1.0 - P;
  97.          B       := A + Y + 1.0;
  98.          Term    := 0.0;
  99.          Pn[ 1 ] := 1.0;
  100.          Pn[ 2 ] := Y;
  101.          Pn[ 3 ] := Y + 1.0;
  102.          Pn[ 4 ] := Y * B;
  103.          Gin     := Pn[ 3 ] / Pn[ 4 ];
  104.          Done    := FALSE;
  105.          Iter    := 0;
  106.  
  107.          REPEAT
  108.  
  109.             Iter   := Iter + 1;
  110.             A      := A + 1.0;
  111.             B      := B + 2.0;
  112.             Term   := Term + 1.0;
  113.             An     := A * Term;
  114.  
  115.             Pn[ 5 ] := B * Pn[ 3 ] - An * Pn[ 1 ];
  116.             Pn[ 6 ] := B * Pn[ 4 ] - An * Pn[ 2 ];
  117.  
  118.             IF( Pn[ 6 ] <> 0.0 ) THEN
  119.                BEGIN
  120.  
  121.                   Rn     := Pn[ 5 ] / Pn[ 6 ];
  122.                   Dif    := ABS( Gin - Rn );
  123.  
  124.                   IF( Dif <= Eps ) THEN
  125.                      IF( Dif <= ( Eps * Rn ) ) THEN Done := TRUE;
  126.  
  127.                   Gin    := Rn;
  128.  
  129.                END;
  130.  
  131.             Pn[ 1 ] := Pn[ 3 ];
  132.             Pn[ 2 ] := Pn[ 4 ];
  133.             Pn[ 3 ] := Pn[ 5 ];
  134.             Pn[ 4 ] := Pn[ 6 ];
  135.  
  136.             IF( ABS( Pn[ 5 ] ) >= Oflo ) THEN
  137.                BEGIN
  138.                   Pn[ 1 ] := Pn[ 1 ] / Oflo;
  139.                   Pn[ 2 ] := Pn[ 2 ] / Oflo;
  140.                   Pn[ 3 ] := Pn[ 3 ] / Oflo;
  141.                   Pn[ 4 ] := Pn[ 4 ] / Oflo;
  142.                END;
  143.  
  144.          UNTIL ( Iter > MaxIter ) OR Done;
  145.  
  146.          Gin    := 1.0 - ( F * Gin * P );
  147.  
  148.          GammaIn := Gin;
  149.                                    (* Calculate precision of result *)
  150.  
  151.          IF Dif <> 0.0 THEN
  152.             Cprec := -LogTen( Dif )
  153.          ELSE
  154.             Cprec := MaxPrec;
  155.  
  156.       END   (* Continued Fraction *)
  157.  
  158.    ELSE
  159.  
  160.       BEGIN (* Infinite series *)
  161.  
  162.          Iter    := 0;
  163.  
  164.          Term    := 1.0;
  165.          C       := 1.0;
  166.          A       := P;
  167.          Done    := FALSE;
  168.  
  169.          REPEAT
  170.  
  171.             A       := A + 1.0;
  172.             Term    := Term * Y / A;
  173.             C       := C + Term;
  174.  
  175.             Iter    := Iter + 1;
  176.  
  177.          UNTIL ( Iter > MaxIter ) OR ( ( Term / C ) <= Eps );
  178.  
  179.          GammaIn := C * F;
  180.                                    (* Calculate precision of result *)
  181.  
  182.          Cprec := -LogTen( Term / C );
  183.  
  184.       END   (* Infinite series *);
  185.  
  186. 9000: (* Error exit *)
  187.  
  188. END    (* GammaIn *);
  189.